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ABSTRACT 



Previous experimental investigations have shown that the characteristics of flow 
about a circular cylinder immersed in a time-dependent flow exhibit cvcle-to-cycle 
variations. These variations have been attributed to the variations in the spanwise 
coherence, aspect ratio, nonuniformity of the flow, and random disturbances in the 
ambient flow. A theoretical investigation was undertaken to examine the stability of 
the How characteristics in terms of the initial state of the vortices. An idealized model 
has been devised and the position of the vortex was varied systematically. The results 
have shown that finite-precision information about the characteristics of the flow does 
not lead to finite-precision information at a later stage. In fact the advection of the 
vortices can give rise to chaotic behavior in the calculated lift and drag forces and in 
the velocity field. It is concluded that the cycle-to-cycle variations are not entirely due 
to lack of spanwise coherence and that they are mostly a consequence of the chaotic 
motion which can result from the advection of the vortices in a time-dependent flow. 
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I. INTRODUCTION 



In the past ten years there has been much effort to analyze the separated time 
dependent flow about bluff bodies in general and circular cylinders in particular. 
Shortfalls in the M orison Equation to adequately predict in-line forces as well as the 
advent of sophisticated computational and experimental techniques have stimulated a 
great deal . of research in this area. The accumulation of this numerical, experimental 
and visual data has generated the need for a deeper understanding of the motion of 
vortices in time dependent flow. Specifically in oscillating flow the stability of the 
vortices generated is of particular importance. The significance of the stability of the 
vortices is a result of its possible exposition of documented cycle-to-cycle variations in 
measured forces on the cylinder. 

In this application, stability is not a characteristic of the vortex life or demise but 
of the vortex motion or path. A system is referred to as stable in any particular state 
of equilibrium if, after being given a finite perturbation, it tends to return to the state 
of equilibrium existing before it was disturbed. On the contrary an unstable system 
subjected to a finite perturbation catastrophically fails to return to its initial state of 
equilibrium. Many real systems in nature do not fall into either of these two categories 
though. In these systems there are certain times when finite disturbances do not lead 
to finite results nor catastrophic or runaway reactions. These complicated systems are 
described as neither stable nor unstable but as systems susceptible to chaos. 

Whether or not the stability of the vortex motion in oscillating flows can be 
categorized as chaotic requires not only an analysis of the vortex motion but also a 
clear definition and understanding of chaos. Unfortunately the pursuit of this 
understanding remains the subject of ongoing research and a concise definition is yet 
elusive. Strides have been made however in identifying chaotic behavior. Systems 
studied that demonstrate chaotic behavior could best be described as bifurcating 
systems that operate along a razor's edge or precipice. The slightest disturbance will 
push the system one way or the other and will lead to decidedly different (but not 
unstable) results. Additionally, the systems studied reveal that the state of chaotic 
behavior is not all encompassing but exists more in regions of susceptibility to chaotic 
behavior preceded and followed by regions of stable behavior. 
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The application of this analysis of chaotic behavior to vortex motion in 
oscillating flow is particularly useful due to the complexity of the factors driving the 
vortex trajectory. The coupling of relative position to other vortices (real or image) in 
-a time dependent flow is both intricate and inseparable. Additionally, the flow field 
about a bluff body near the surface, combined with the problems of multiple vortices, 
viscous dissipation, wall reflection and boundary layer and limited coherence length, is 
extremely intricate. Consequently, the motion of the vortex at times becomes subject 
to finite perturbations without finite results. The specific coupling that makes the 
system susceptible to chaotic behavior gives rise to the regions of chaos preceded and 
followed by regions of stability. 

In order to understand these complicated causes and effects first an 
understanding of the motion of a single potential vortex in oscillating flow about a 
bluff body is required. Therein lies the impetus of this study. In an ideal two 
dimensional potential flow, wall reflection, wall boundary layers and coherence length 
are eliminated. Additional vortices can be added to study mutual induction effects. 
Further, advanced models of the vortex (i.e. Rankine and Gaussian time dependent 
models) can clarify the broad effects of viscous dissipation. 

The scope of this study includes the potential flow of one and two vortices placed 
in the vicinity of a circular cylinder in oscillating flow. Emphasis is placed on 
identifying the onset of chaotic behavior and its effects on the flow characteristics. 
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II. BACKGROUND 



The concept of chaos in dynamic systems is by no means new. However, the 
detailed investigation as to the onset and mechanisms of chaos is still comparatively 
young. The primary difficulty in the investigation of chaos is the fact that by definition 
chaotic behavior is non-integrable. Because of this the first attempts to define the 
concept of chaos were undertaken by mathematicians. In 1984 Holmes [Ref. lj 
investigated chaotic motion in forced oscillations. Here the onset of chaos was linked 
to periodicity or aperiodic behavior of the function through scrutiny of the harmonics. 
Although this is a credible definition its application appears limited to the idealized 
mathematical system from which it was developed. In 1980 Aref [Ref. 2] expanded on 
this definition and the use of Poincare' maps by illustrating chaotic behavior in paths 
of finite, ideal, and controllable vortices. His motivation for the investigation was the 
onset and mechanisms of chaos. Therefore, while the idealized system that he designed 
was useful for demonstrating chaotic behavior, its application to real world dynamics 
was severly limited. 

Other efforts were made to try to link the study of chaos to systems more closely 
resembling those found in nature. In 1984 Hardin & Mason [Ref. 3] investigated the 
two dimensional system of axial vortices in a pipe. This was the first attempt to 
examine a system that resembled real phenomena and that also underwent chaotic 
behavior. All previous studies examined chaotic behavior by generating systems 
susceptible to chaos and then extrapolating the results to possible real world problems. 

Finally, in 1986 Dowell & Pezeshki [Ref. 4] investigated the one dimensional 
Duffings Equation and the onset of chaos. In its relation to real systems Duffings 
Equations could be used to describe the snap buckling of a beam and other 
phenomena. In the beam the onset of chaos is easily observed but the integration of 
the equation to predict or study the mechanisms of chaos is somewhat abstract. Of 
particular interest in this study was the examination of the onset of chaos by Phase 
Plane Plots (displacement - velocity plots). 

By reasons of their treatment of fluid dynamics in general and vortex motion in 
particular, Arefs chaotic advection and Hardin & Mason's vortices in a pipe have 
direct application to this investigation and will be discussed in greater detail. 
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A. CHAOTIC ADVECTION 

Aref [Ref. 2] devised a numerical model to represent the motion of an 
incompressible inviscid fluid within a circular bounding contour. A stirring agitator 
was modelled as a potential point vortex. The agitator(s) provided the potential flow 
within the boundary. A global Eulerian view of an array of markers was traced to 
examine and verify the poor efficiency of the laminar mixing. 

In a follow-up study, he examined the case of two fixed agitators (point vortices) 
in the fluid that generated "piecewise-constant stirrer motion." The vorticity of the 
agitators was controllable (on or off) and switched back and forth at prescribed times 
to generate unsteady flow. Again a marker array of particles was traced. The mixing 
was clearly superior and Aref classified the motion as chaotic. The parameters that 
dictate the switching of the vortices and therefore the chaotic regions were cited as 
mechanisms for chaos. 

Extrapolations to real phenomena include enhanced mixing chambers with 
controllable agitators and similarities to gravitational systems. 

B. VORTICES IN A PIPE 

The motion of axial vortices was investigated by Hardin and Mason [Ref. 3] 
because of its widespread appearance in physical situations. These vortices occur in 
virtually every’ flow that changes axial directions (i.e. rivers, pipes,. lung bronchial tubes, 
etc.). In these physical systems pairs of counterrotating vortices are produced. This 
study investigated the effects of one and two pairs of equal strength counterrotating 
two dimensional potential vortices enclosed in a circular boundary. Vortex advection 
by induced velocities (real or image) was determined and displayed by means of 
Eulerian paths of the vortices. 

Controlling factors on this investigation were initial locations of the vortices 
within the boundary'. For the case of single pair of vortices the paths were found to be 
periodic or stable. In the case of two-vortex pairs (four vortices within the boundary), 
the periodicity of the vortex path and flow was a function of where the vortices were 
placed initially. This supports the concept of chaotic motions occurring within regions 
of susceptibility. It w r as found that within these regions of susceptibility that small 
perturbations (0.001 radii of the enclosing boundary) were sufficient to drive the vortex 
path aperiodic and demonstrate chaotic behavior. Hardin & Mason [Ref. 3] do note 
that experimental evidence indicates that bifurcating ■ or branching flows are 
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quantitatively more stable in real viscous flows than these numerical results suggest. 
This is attributed to the complexity of this secondary flow as well as viscous 
dissipation, vortex decay and coring influences. It does not, however, diminish the 
validity of the study or the qualitative results of the susceptibility of fluid flows to the 
onset of chaos. 
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III. ANALYSIS 



A. INTRODUCTION 

In 1976 Sarpkaya [Ref. 5] presented extensive data on oscillating flow about 
circular cylinders. In addition to modifying the long standing M orison Equation, this 
data also provided extreme insight into transverse forces on the cylinder. Performed at 
Reynolds numbers of prototype magnitude and across a wide spectrum of Keulegan- 
Carpenter [Ref. 6] numbers, these experiments provided a data base to which all 
subsequent publications and literature on the subject has been based. Of particular 
significance and not predicted by the M orison Equation was the documentation of 
cycle-to-cycle variations in pressures and forces across the cylinder in pure harmonic 
flow (see Figure E.l). These variations, particularly significant in the transverse 
direction, were attributed primarily to limited coherence lengths and to a lesser extent 
wall reflections and boundary layers. 

The distorting effects of coherence length on the wake of a circular cylinder in 
uniform flow was clearly discussed by Gerlach [Ref. 7] and Graham [Ref. 8.] In 
oscillating flows, however, the von Karman street forms perpendicular or transverse to 
the flow direction in the range 8<K < 13. This makes experimentation to identify the 
coherence length extremely difficult. In an effort to specifically examine the effects of 
spanwise coherence, O'Keefe [Ref. 9] in 1986 undertook similar experiments to 
Sarpkaya [Ref. 5] but with the minimum aspect ratio (L/D = 2.0) that could be 
accurately instrumented. These nearly two dimensional experiments provided 
extremely accurate data that was not blurred by the influence of a limited coherence 
length. Cycle-to-cycle variations in in-line and transverse forces remained though; 
particularly in the Keulegan-Carpenter number range 8 < K c < 13. 

Knowing that at any given point in time the forces acting on the cylinder is 
governed primarily by one dominant vortex, a detailed investigation of one potential 
vortex in the vicinity of a circular cylinder in harmonically oscillating flow was dictated. 
After an investigation of the single vortex case the follow-up was to investigate two 
potential vortices in the vicinity of a circular cylinder. 
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B. POTENTIAL FLOW MODELING 
1. One Vortex Model 

In the complex plane the potential How field for a vortex in the vicinity of a 
circular cylinder in harmonically oscillating flow is given by the following expression: 

fa) = U max sin((0t X z + z 1 ) + i f" ln ( z_z o) -i f’ ln ( z_ = _ ) (3-0 

27t ° 2n z 

o 

where: c - radius of the cylinder 

T - strength of the vortex 
CO - frequency of the oscillating flow 
U max - maximum amplitude of the flow velocity 
z Q - location of the vortex 

In nondimensional form this expression becomes: 

f (0 = sin(27ttXC + -}-) + iKln(C ~ t 0 ) - iKlntf - (3.2) 

where: K.= r/(27tcU max ) 

t = t/T (T - period of the flow) 

C, = z/c 

This expression can be differentiated with respect to C, to determine the flow 
field velocities. At q = q o the expression for the velocity of the vortex becomes: 

Ifl - u- iv - sin(27tt)(l -4 t) - (3.3) 

S s=; o v 

where r is the radius between the real and image vortex and is given bv: 
v ® 




(3.4) 



(Note that the second term of eqn 3.2 is dropped since a vortex does not 
impart velocity to itself.) 

Calculations for the forces acting on the cylinder were made by integrating the 
Blasius Force Equation: 



D-iL = — ip$(— ) 2 dz + ip— If dz (3.5) 

2 oz dt 

Nondimensionalized and evaluated at ^ = this equation reduces to the 
following form: 

D'-iL' - ^-cos(27tT) + £-(C-±.) (3.6) 

k k oz 0 r 

C C tjO 

where: K c = (U max T)/(2c) 

2. Two Vortex Model 

Invoking the same procedure as the one vortex model it can be shown that the 
velocity of two vortices A and B in the vicinity of a the circular cylinder is given by: 



QL\ = u-iv 






= sin( 2 jtt)( I -■ 



A 



72 )_( 



iK 



Sa“ 1 /? 



—)-{■ 



iK 



Ca-Sb 



) + ( 



iK 



5 A - ^ 



-) (3-7) 



B 



— | = u — iv 



sin( 2irt)( I — 4) - 1 ; — iiy-) - (— & — ) + !; — 

ft - '/r A ft “ft ^D _ 'A B 



■) (3.8) 



3. Rankine Vortex Model 

A standard rankine vortex model was used with linear velocity distribution 
between 0 and r core . The value of the core radius was chosed to be 0.1(1 times the 
radius of the cylinder generating the vortices. 
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4. Gaussian Time Dependent Vortex Model 

A standard Gaussian time dependent vortex model was used where the 
generalized expression for the velocity distribution is: 



V = — 2 (1 ~ exp( 



27tr 




(3.9) 



Here t is not the clock time of the flow but the life time of the vortex. It has 

v 

been shown that differentiating this with respect to the radius yields r = 2.24V(vt v ) as 
the expression for the core radius of the vortex. Nondimensionalizing eqn 3.8 in terms 
of the oscillating flow parameters, and examing the velocity of a vortex from eqn 3.3. 
the velocity of a Gaussian time dependent vortex is given by: 



For the modeling the Reynolds Number was chosen to be Re= 10,000 and t 
was set such that the core of the vortex was r =0.10 at inception. 

5. Numerical Integration and Accuracy 

In this analysis a finite differencing technique was employed. The study of 
potentially chaotic systems requires a high degree of accuracy. Due to its guaranteed 
stability and overall efficiency, the Adams-Moulton predictor corrector method was 
selected to convect the vortices: 



u — iv = sin(27CT)( 1 — — ^ 



1 \ _ iK 



{1 - exp( - 



r v V(K c /Re ) 
16(t v /T) 



•)) (3-10) 



-np 



= z 



n-1 



+ (55 q n .j - 59 q n _2 + 37 q n . 3 - 9 q n _ 4 ) 



At 
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(3.11) 



z n = z n-l + < 9 9 np + 19 fln-l " 5 %-2 + 9n- 3 ) < 3 - 12 ) 

In addition to ensured stability, this method provided fourth order accuracy 
with local error 0(h 5 ). With 1440 time steps per cycle, the local error remained at the 
limits of double precission machine accuracy. Initial startup was made using the fourth 
order Runge-Kutta scheme. All continuations were based on the last four calculated 
locations of the vortices, thereby ensuring accuracy and duplication at any cycle. 
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While discussing accuracy it must be noted, as mentioned earlier, that by 
definition chaotic behavior is non-integrable. Whatever scheme used can only illustrate 
qualitatively the onset of chaos. If used for comparison universally, all methods yield 
equally valid qualitative results. The effort devoted here to accuracy is therefore to 
ensure the validity of the results and add to overall integrity of this report. 
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IV. RESULTS 



A. GENERAL 

From Sarpkaya [Ref. 5] and other studies, it is concluded that the region of 
special interest in the field of oscillating flows is confined to 8 < K c < 13, an interval in 
which a half transverse vortex street is generated. It is because of this reason that in 
the present numerical experimentation K c was chosen equal to 10.0. 

The vortex path is a function of both the imposed ambient flow field and the 
velocities induced by the vortices (real or image). This coupling was first investigated 
through the use of a single vortex in uniform flow in the vicinity of a circular cylinder. 
Figure E.2 illustrates the path of the vortex for various vortex strengths (K = 0.5, 0.625, 
0.75). It is seen that the capture by the image vortex (partial looping) lasts only until 
the flow field around the cylinder pushes the vortex away from the image. It is also 
seen that this capture is a function of the initial location of the vortex. When the 
vortex is captured, sharp cusps appear where the direction of motion of the vortex is 
quickly reversed. These characteristics played even more decisive roles in the unsteady 
case. With the exception of the foregoing, the vortex strengths used in the numerical 
experimentation were all set at K = 0.50. 

B. ONE VORTEX MODEL 

Initial experimentation included varying the starting locations and times of the 
single point vortex in order to observe the vortex path and force coefficients. Figures 
E.3 through E.6 show the force coefficients for various starting locations and times and 
varying numbers of cycles. These figures clearly show cycle-to-cycle variations in both 
the inline and transverse force coefficients. In comparison to experimental force traces 
shown in Figure E.l, of comparable K c , the variations are of similar magnitude and 
shape. 

Understanding the origin of these cycle-to-cycle variations required a scrutiny of 
the vortex paths themselves. Figures E.7 and E.8 show two starting positions of the 
vortex of (0, 1.5083) and (0, 1.5084). The oscillating flow velocity was set at L' o = 0 + 
(superscript indicates the initial x-direction). Figures E.7(a) and E.8(a) show that for at 
least the first five cycles of the flow the vortex paths are virtually identical. Also noted 
in this figure are sharp cusps where the overall direction of the vortices is abruptly 
changed. 
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Although the difference in the begining locations of the sixth cycle is virtually 
inperceptible = (0.995, 0.723) vice ^ = (0.970, 0.753)}, the vortex in Figure E.7(b) 
forms a sharp cusp at the top of the cylinder (geographical location of the maximum 
velocity at any time in the flow) but the vortex in Figure E.8(b) is captured by its 
image and orbits the cylinder. The ending points for cycle six are grossly different (see 
Figures E.7(c) and E.8(c)} and the vortex paths never again become similar. The finite 
perturbation in starting location has yielded explicitly different results. Nevertheless, 
the vortex paths in Figures E.7(c) and E.S(c) remain bounded. Extending the 
calculations to 100 cycles of the flow, Figure E.9 shows that while neither vortex paths 
are identical, both remain bounded. They clearly seem to slide in and out of chaotic 
behavior, preceded and followed by periods of stability. 

Figure E.10 is a plot of the magnitude of the velocity of the vortices. It is again 
evident that these two vortices have virtually identical velocity magnitudes until the 
sixth cycle of flow. The capture and looping depicted in Figures E.7(b) and E.S(b) is 
also evident in Figure E.10. Breaking down this velocity magnitude into its u and v 
components, Figures E. 11 through E.13 show that overall periods of stability are 
intermixed with these vortex captures. 

Figures E.14 through E.19 and E.20 through E.25 are similar investigations that 
document two starting points of the vortices where both initial x-locations were set at 
x=0.0 and the oscillating flow velocity was set at U 0 = U+ ax and U o = U“ ax , 
respectively. Figures E.26 through E.31 are similar investigations that document two 
starting points of the vortex in the vicinity of x= 1.0, y= 1.0, and the oscillating flow 
velocity was set at U o = 0 + . These investigations support the conclusions of the first 
investigation and demonstrate once again a finite precision input does not yield a 
finite precission output. They also demonstrate the periods of stability that precede 
and follow chaotic behavior, (i.e., chaos in order and order in chaos). 

As seen in the proceeding set of experiments, the looping or continued looping of 
the vortex tightly around the cylinder, was a function not only of the location of the 
vortex but also of the current state of the oscillating flow. This suggests that at a 
given initial location, a variance of initial velocity or phase shift of the oscillating flow 
could drive the system from stable to chaotic behavior at a given point in time. Figure 
E.32 is a series of vortex paths that supports this concept. Four vortex paths were 
generated with the initial location of x = y= 1.0 and a variance of initial velocity by a 
phase shift from 0.75 (U 0 =U}^ ax ) to 1.0 (U o = 0 + ). It is clear that while the Figure 
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E. 32(a) documents stability even after ten cycles, the cusps and asymmetric tracks in 
Figure E. 32(d) demonstrate the mixture of stable and chaotic regions. 

Figure E.33 depicts the results of a similar investigation where the initial location 
of the vortices is at (0, 2.2192) and (0, 2.2193). These results also show that the 
chaotic behavior depends not only on the vortex location but also on the flow 
environment. In fact, the two are coupled inseparably. 

The stability and symmetry of Figures E. 32(a) and E. 33(a) also suggest, however, 
that chaotic motion may not be inevitable. They infer that the location of the vortex 
and field flow may be tuned in such a way as to permanently avoid chaos. This would 
be paramount to permanently balancing on a razor's edge. However, such a balance is 
unstable. As seen in the first investigation, the slightest disturbance is sufficient to 
push the -system into chaotic behavior. 

The final phase of the one vortex experiment was to examine the effects of higher 
order models of a potential vortex in an effort to better understand vortices in nature. 
Figure E.34 is a comparison of the results obtained with ideal, Rankine, and Gaussian 
time dependent vortex models. It is clear from these figures that while the onset of 
chaos may be different, the chaotic behavior and susceptibility to chaos remains 
virtually the same. 

C. TWO VORTEX MODEL 

Separated flow past a bluff body generates pairs of counterrotating vortices. In 

both steady and oscillating flow about a circular cylinder, the first vortex pair is formed 

symmetrically and is significantly stronger than subsequent asymmetric vortices. 

Additionally, it has been documented that in oscillating flow in the region 8<K c < 13, 

a half von-Karman vortex street forms perpendicular or transverse to the flow 

direction. In this section of the experimentation, two vortices -were placed in the 

vicinity of a circular cylinder immersed in a harmonically oscillating flow. The vortices 

were tracked through various numbers of cycles of flow. The locations and strengths 

of the vortices were varied in order to examine the mechanism for the formation of a 

half von-Karman street. Additionally, the stability of the vortices as discussed in the 

one vortex model was examined. For all of the investigations, the characteristic of the 

flow was set at K c = 10.0. The strength of the vortices was set at K = 0.50±5 

(although contrary in direction), and their locations at C^ = (a,a) and 

Ug = (a±5,-a±5), where 0.80<a < 1.414, and 5< < 1. Based on flow visualization, the 

initial velocity was set at U = U“ Y . 

' o max 
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The first part of the investigation required the analysis of vortices, symmetric in 
strength and location. Figure E.35 demonstrates that within the regions examined, the 
mutual induction velocities between the real vortices dominate the advection of the 
vortices. It is also clear that the vortices are quickly convected away in parallel paths. 
The initial radial location determines the separation distance of the two parallel paths. 
This separation distance combined with the strength of the vortices determines the rate 
at which the vortices are swept away. 

Figures E.36 through E.38 are of a series of investigations where the vortices 
have equal strengths (K^ = Kg) but the locations are asymmetric (|^g| = £ r x 
The overall result is a general curving of the vortex path toward the vortex that had 
the larger initial radial distance from the cylinder. This is a result of the image vortices 
having a greater effect on the vortex that was closer to the cylinder. It is seen from 
Figure E.37 that a 10% difference (c r = 0.90) in the two starting locations is sufficient 
to significantly curve the vortex paths. The vortex paths do not gravitate to each other 
though and the vortices do not become coincident. 

Figures E.39 through E.42 are of a series of investigations where the vortices 
have symmetric initial locations (^ = Cg) anc * asymmetric strengths (Kg = £ K ><K^). 
The result here is a curvature similar to that of the previous investigation. From 
Figure E.41 it is seen that only a 5% difference (C K = 0.95) in the vortex strengths is 
sufficient to significantly alter the vortex paths. Clearly the results indicate that the 
system is much more sensitive to asymmetry in vortex strengths than to asymmetry in 
vortex location. Again the vortices do not become coincident. 

Figure E.43 is a combination of the previous investigations where there is 
asymmetry in both the vortex strengths (Kg = e K ><K^) and in initial location 
(|^g| = £ r x I^aD- Again both vortex tracks bend toward the dominant vortex but as 
expected, with a more severe curvature. Although in real flow two vortices are 
generated every time the oscillating flow direction changes, this study of two vortices 
moving together to either the top or bottom of the cylinder is particularly significant. 
It suggests the evolution of a half von-Karman vortex street transverse to the flow 
based on the asymmetry of the vortices in the street. 

Finally, an investigation was made where the asymmetry of the vortices was so 
extensive that the vortex paths curved back to the vicinity of the cylinder. Although in 
viscous flow the life of a vortex is finite, combinations of additional vortices coupled 
with the flow field make this situation plausible and warrants its investigation. Figure 
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E.44 documents the sharp cusps, capture, and chaotic behavior of the vortices as 
discussed in the previous section. With the additional vortex (and image) the coupling 
is much more complicated between the geographical location of the vortex and the 
flow field. The identification of chaotic behavior is none the less possible. 
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V. CONCLUSIONS 



The results presented herein warranted the following conclusions: 

1. Finite precision information about the initial characteristics of a vortex (or of 
larger number of vortices) provides no finite-precision information about the 
characteristics of flow at later times and that the lack of spanwise coherence of 
the vortices is not the sole cause of cycle-to-cycle variations in the measured 
quantities. 

2. An idealized vortex model suggests that the sensitivity of the characteristics of 
the flow to random disturbances superimposed on the motion of vortices at an 
earlier time is the main flow feature responsible for the observed cycle-to-cycle 
variations. 

3. The chaotic behavior of the vortex in oscillating flow occurs in regions of 
susceptibility similar to those found by Dowell and Pezeshki [Ref. 4] with 
similar sensitivity. 

4. The susceptible regions for chaos in oscillating flow are coupled inseperably 
with the location of the vortices and the state of the environmental flow field. 
In other worlds, for a given flow phase there is a vortex location that will lead 
to chaos. Likewise, for every vortex position near the cylinder, there is a 
particular flow phase that will eventually lead to chaos. 

5 . The formation of a half von-Karman vortex street in the range of 8<K c < 13 

may be attributed to the asymmetry in the formation of the vortices. 

6. The extent to which the sensitivity of the flow to the conditions prevailing in 

the previous cycles is damped by finite core effects and viscosity has been 

investigated and it has been shown that they change only the time of onset of 
the chaos, not its invevitable occurrence. 
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APPENDIX A 

ONE VORTEX COMPUTER MODEL 



c ********************************************************************* 



c * * 

c * VORTEX * 

C * * 

C * BY LT. W.T. MCCOY * 

C * * 



Q ************************************************************^^^:^r^^:^:^:^: 

c 

C DESCRIPTION 

C 

C THIS PROGRAM TRACKS A VORTEX IN COMPLEX SPACE AROUND A UNIT CYLINDER 
C IN OSCILLATING FLOW. INITIAL LOCATION, STARTING TIME, AND THE NUMBER 
C OF CYCLES OF THE OSCILLATING FLOW ARE REQUIRED INPUT. 

C 

C THE PROGRAM IS DESIGNED TO RUN ON DISSPLA BY USING THE FOLLOWING 
C DISSPLA EXEC COMMAND : 'DISSPLA VORTEX'. IT IS CURRENTLY 
C DIMENSIONED TO RUN A MAX OF 10 CYCLES AT 1440 STEPS PER CYCLE. 

C 2.0M OF VIRTUAL MEMORY IS REQUIRED TO RUN THE PROGRAM AT MAXIMUM 
C CAPACITY. THE PROGRAM IS INTERACTIVE FOR INPUT OPTIONS. 

C 

C ALL CONSTANTS ARE SET IN THE MAIN PROGRAM. PREVIOUS INPUT OPTIONS 
C ARE READ FROM FILE 'FILE FT25F001 ' . FAILURE TO GENERATE THIS FILE 
C PRIOR TO INITIAL RUNNING RESULTS IN A ERROR MESSAGE AND TERMINATION 
C OF THE RUN. PRIOR TO RUNNING, FILE 'FILE FT25F001' SHOULD BE 
C INITIATED CONTAINING 4 COMPLEX NUMBERS , 1 REAL NUMBER, AND 1 INTEGER. 
C WITH THIS INPUT THE PROGRAM DOES NOT NEED TO BE TO BE RECOMPILED 
C FOR EACH RUM. 

C 

C THIS PROGRAM IS WRITTEN IN DOUBLE PRECISION AND EMPLOYS THE MILNE 
C PREDICTOR-CORRECTOR CONVECTION SCHEME, WITH INITIAL POINTS GENER- 
C ATED BY A 4TH ORDER RUNGE-KUTTA SCHEME. 
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c 

C OUTPUT OF THE PROGRAM INCLUDES THE LAST FOUR VALUES OF THE LOCATION 
C OF THE VORTEX, FINISH TIME, AND THE NUMBER OF CYCLES. THESE VALUES 
C ARE WRITTEN IN FILE 'FILE FT25F001 1 AND ALLOW FOR OBSERVATIONS BEYOND 
C THE MAX LIMIT OF 10 CYCLES BY CONTINUING WITH ADDITIONAL RUNS. 

C 

C OUTPUT FOR THE VORTEX PATH AS WELL AS FORCE COEFFICIENTS AND 
C VELOCITY MUNI PULAT IONS ARE MADE GRAPHICALLY. SELECTIONS OF PLOTS 
C FOR OUTPUT ARE MADE INTERACTIVELY. PLOTS MAY BE MADE ON 
C A TEK618 OR SENT TO A METAFILE (PRE-EDIT SUBROUTINE 'GRAPH'). 

C 

C 

c 

c 

C ALPHABETICAL LISTING OF VARIABLES 

C 

C VARIABLE (SUBROUTINE USED) DESCRIPTION. 

C 

C 

c 

C CD(N) (GRAPH) NON-DIMENSIONAL ARRAY OF INLINE (DRAG) 

C FORCE COEFFICIENTS FOR VORTEX POSITION Z(N) . 

C 

C CL(N) (GRAPH) NON- DIMENSIONAL ARRAY OF PERPENDICTULAR 

C (LIFT) FORCE COEFFICIENTS FOR VORTEX POSITION 

C Z (N) . 

C 

C COEF (FORCE, GRAPH) COMPLEX VARIABLE CONTAINING AS ITS 

C PARTS THE NON-DIMENSIONAL FORCE COEFFICIENTS. 

C 

C DELTAT (STEP, Q, OUTPUT, FORCE, GRAPH) INCREMENTAL STEP 

C IN NON-DIMENSIONAL TIME (INVERSE OF STEPS PER 

C CYCLE). 

C 

C DELZ (STEP ) COMPLEX INCREMENTAL STEP OF THE LOCATION 
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C OF THE VORTEX (FROM Z1 TO Z2) . 

C 

C DFDZ (Q)COMPLEX PARTIAL DERIVATIVE OF THE COMPLEX 

C FUNCTION F WITH TO RESPECT Z. 

C 

C DISP(N) (STEP, GRAPH) ARRAY CONTAINING THE TOTAL DISPLACEMENT 

C OF THE VORTEX FOR A CORRESPONDING LOCATION Z(N). 

C 

C I (Q, FORCE) CONSTANT COMPLEX NUMBER = (0.0, 1.0) 

C 

C ILHS (GRAPH) LEFT HAND VALUE OF THE TIME AXIS OF THE 

C FORCE COEFFICIENT PLOT. 

C 

C IRHS (GRAPH) RIGHT HAND VALUE OF THE TIME AXIS OF THE 

C FORCE COEFFICIENT PLOT. 

C 

C ISTEPS... (ALL) NUMBER OF INCREMENTAL STEPS PER CYCLE TO BE 

C TAKEN FOR ALL ITERATIONS. 

C 

C Kl-4 (RUNGE) RUNGE-KUTTA CONSTANTS 

C 

C KAPPA (ALL) STRENGTH OF THE VORTEX (GAMMA/ (2*PI*UMAX)) 

C 

C KC (ALL) KEULLEGAN- CARPENTER NUMBER FOR THE OSCILLATING 

C FLOW. 

C 

C LI (INPUT, GRAPH) LOGICAL OPERATOR. 

C 

C MUM CYC ( INPUT) TOTAL NUMBER OF CYCLES TO BE EVALUATED. 

C 

C NUMPTS (ALL)TOTAL NUMBER OF POINTS TO BE DETERMINED. 

C (NUMCYC * ISTEPS) 

C 

C PHSHFT (GRPAH) INITIAL PHASE SHIFT OF THE FLOW (0<PHSHFT<1) 

C 
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C PI (ALL) CONSTANT = 3.141592654 

C 

C Q (Q) FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 

C VELOCITY OF THE VORTEX. 

C 

C QMAG(N) (GRAPH) ARRAY CONTAINING THE MAGNITUDE OF THE 

C COMPLEX VELOCITY QVECT(N) . 

C 

C QVECT(N) (STEP, GRAPH) COMPLEX ARRAY CONTAINING THE VELOCITY 

C OF THE VORTEX AT ANY POINT Z(N). 

C 

C RV (Q) VECTOR BETWEEN THE VORTEX AND ITS IMAGE 

C WITHIN THE CYLINDER 

C 

C T (GRAPH) NON- DIMENSIONAL TIME (0-1 FOR A CYCLE). 

C 

C TE (RUNGE,Q) NON- DIMENSIONAL ENVIRONMENTAL FLOW TIME. 

C 

C TT ( STEP, Q,FORCE)NON-DIMENSIONAL TIME CHANGED TO 

C RADIANS FOR COS/SIN EVALUATIONS (0-2*PI FOR A CYCLE). 

C 

C TTF (OUTPUT) NON- DIMENSIONAL FINISH TIME OF THE RUNS. 

C 

C TTI (ALL) INITIAL NON-DIMENSIONAL STARTING TIME FOR THE 

C OSCILLATING FLOW. 

C 

C TV (Q) NON-DIMENSIONAL TIME LIFE OF THE VORTEX. 

C (TE + CONSTANT FOR R-CORE AT INCEPTION) 

C 

C U ( GRAPH) X-DIRECTION VELOCITY MAGNITUDE. 

C 

C V (GRAPH) Y-DIRECTION VELOCITY MAGNITUDE. 

C 

C Z (ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 

C 
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C ZI (GRAPH) IMAGINARY (Y) LOCATIONOF THE VORTEX. 

C 

C ZP (Q) COMPLEX PREDICTED VALUE FOR Z(N) 

C 

C ZR (GRAPH) REAL (X) LOCATION OF THE VORTEX. 

C 

C ZZ ( Q ) DUMMY ARGUMENT FOR THE 'FUNCTION SUBROUTINE Q 

C WHICH CONTAINS THE PASSED VALUE OF Z(N). 

C (RUNGE) DUMMY ARGUMENT CALL FUNCTION SUBROUITNE Q 

C 

c 

c 

c 



C * MAIN PROGRAM * 

c 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPA 

REAL*4 DISP( 14409) ,CD(14409) ,CL(14409) 

COMPLEX*16 Z ( 14409 ) ,QVECT( 14409) . 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC, KAPPA 
1 /FOURTH/ PI 

C 

C SET CONSTANTS 
C 

KAPPA = 0.5D0 
KC = 10. DO 
PI = 3 . 141592654D0 
ISTEPS = 1440 
C 

10 CALL INPUT (Z) 

CALL STEP (Z,QVECT,DISP) 

CALL OUTPUT (Z) 

CALL FORCE (Z, CD, CL) 
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CALL GRAPH (Z , QVECT , DISP , CD , CL) 



C 

WRITE (6 , *) 1 DO YOU WISH TO RE-RUN THE PROGRAM? 1-YES, O-NO' 

READ (5 ,*) LI 
IF (Ll.NE.l) GOTO 20 

WRITE (6,*) 'EXECUTION CONTINUES 1 

GOTO 10 
20 CONTINUE 
C 

CALL DONE PL 
STOP 
END 
C 
C 

Q ********************************************************************* 

C * SUBROUTINE INPUT * 

c ********************************************************************* 
C 

C THIS SUBROUTINE INTERACTIVELY INPUTS THE VALUES FOR THE INITIAL 
C LOCATION OF THE VORTEX, STARTING TIME, AND THE NUMBER OF CYCLES 
C TO BE EXAMINED. 

C 

SUBROUTINE INPUT (Z) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX*16 Z(14409) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 
C 

C READ THE LAST 4 POINTS FOR CONTINUATION OPTION 
C 

REWIND 25 
DO 10 N = 1,4 
10 READ (25 , *) Z(N) 

READ(25 ,*) TTI 
READ (25 , *) NUMCYC 
C 
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C CONTINUATION OPTION 
C 

WRITE (6,*) 'DO YOU WISH TO CONTINUE THIS RUN FORM THE LAST? 1-YES 
# O-NO' 

READ (5 ,*) LI 
IF (L1.NE.0) GOTO 20 

WRITE (6,*) 'INPUT THE NEW STARTING POINT Z(l): (REAL , IMAG) 1 
READ (5 ,*) Z(l) 

WRITE (6, *)' INPUT THE NEW VALUE OF TTI : 1 
READ (5 , *) TTI 
C 

C CALL SUBROUTINE RUNGE TO GENERATE THE FIRST 4 POINTS 
C 

CALL RUNGE (Z) 

C 

C INPUT THE NUMBER OF CYCLES TO BE ANALYZED 
C 

20 WRITE (6,*) 'THE LAST VALUE FOR NUMCYC =',NUMCYC 
WRITE (6 ,*) 1 DO YOU WISH TO CHANGE IT? 1-YES O-NO' 

READ (5 ,*) LI 
IF(Ll.NE.l) GOTO 30 

WRITE (6,*) 'INPUT THE NEW VALUE OF NUMCYC:' 

READ ( 5 , *) NUMCYC 
30 CONTINUE 
C 

NUMPTS = NUMCYC * ISTEPS + 4 
C 

RETURN 

END 

C 

C 

(2 kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k'k'k 

C * SUBROUTINE RUNGE * 

Q kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkikkkkkkkkkkkk'kkkkk'kkk'kkkk 

c 
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C THIS SUBROUTINE USES A 4TH ORDER RUNGE-KUTTA METHOD TO GENERATE THE 
C FIRST FOUR VALUES OF Z. 

C 

SUBROUTINE RUNGE (Z) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX*16 Z (14409 ) , Q ,K1 ,K2 ,K3 ,K4 ,ZP1 , ZZ 
COMMON /FIRST/ TTI,ISTEPS 
C 

DELTAT = 1.D0/DFL0AT(ISTEPS) 

C 

C RUNGE-KUTTA EVALUATION FOR FIRST FOUR POINTS 
C 

DO 5 N = 1,3 
ZZ = Z(N) 

TE1 = (DFLOAT(N-l) * DELTAT + TTI ) 

K1 = DLETAT * Q(ZZ,TE1) 

ZZ = Z(N)+0 . 5D0*K1 
TE2 = TE1+ 0 . 5D0*DELTAT 
K2 = DELTAT * Q(ZZ,TE2) 

ZZ = Z(N)+0.5D0*K2 

K3 = DELTAT * Q(ZZ, TE2) 

ZZ = Z(N) + K3 

TE4 = TE1 + DELTAT 

K4 = DELTAT * Q(ZZ, TE4) 

C 

5 Z(N+1) = Z(N) + (K1-2.D0*K2+2.D0*K3+K4)/6.D0 
C 
C 
C 

RETURN 

END 

C 

C 

c ********************************************************************* 

C * SUBROUTINE STEP * 
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(2 xxxxxxxx:K*xx:**x;kx***xx:**:****;A;*****:**x 

c 

C THIS SUBROUTINE PROGRESSIVELY CALCULATES THE VALUES OF Z(N), BASED ON 
C A PREDICTOR-CORRECTOR CONVECTION SCHEME. 

C IT ALSO DETERMINES THE VELOCITY VECTOR, MAGNITUDE AND PARTS FOR 
C FUTURE MUNI PULAT I ON AS WELL AS OVERALL DISPLACEMENT MAGNITUDES. 

C 

SUBROUTINE STEP (Z,QVECT,DISP) 

IMPLICIT REALMS (A-H,0-Z) 

REAL *4 DISP(14409) 

REAL*3 KC, KAPPA 

COMPLEXES Z(14409) ,QVECT(14409) ,DELZ,ZP ,Q,ZZ 
COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS/THIRD/ KC, KAPPA 
C 

DELTAT = 1 . DO/DFLOAT ( I STEPS ) 

C 

C EVALUATE VELOCITIES AND DISPLACEMENTS FOR FIRST 4 POINTS 
C 

DISP(l) = 0.0 
C 

DO 7 N = 1,4 

TE = (DFLOAT(N-l) * DELTAT + TTI ) 

ZZ = Z (N) 

QVECT(N) = Q(ZZ,TE) 

7 DISP (N+l ) = Z(N+1 ) - Z (N) 

C 

DO 10 N = 5, NUMPTS 

NM1 = N-l 

NM2 = N-2 

NM3 = N-3 

NM4 = N-4 

TE = (DFLOAT (N-l ) *DELTAT+TTI ) 

C 

C CALCULATE PREDICTED VALUE OF Z(N) 

C 
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ZP = Z(NM1)+(55.D0*QVECT(NM1) -59 .D0*QVECT(NM2) +37 .D0*QVECT(NM3 ) 
1 -9.D0*QVECT(NM4) )*DELTAT* 2. DO * KC / 24. DO 
C 

C CALCULATE ACTUAL VALUE OF Z(N) 

C 

DELZ = (9 .DO*Q(ZP,TE) + 19 . D0*QVECT (NM1 ) -5 .D0*QVECT(NM2) 

1 +QVECT(NM3) )*DELTAT * 2. DO * KC / 24. DO 
Z(N) = Z(NH1) + DELZ 
C 

C EVALUATE VELOCITIES AT Z(N) 

C 

QVECT(N) = Q(Z(N) , TE) 

10 DISP(N) = DISP(NMl) + ABS (DELZ) 

C 

C 

RETURN 

END 

C 

C 

C - - 

C | FUNCTION SUBPROGRAM Q | 

C - - - 

c 

C THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF Q 
C - 

COMPLEX FUNCTION Q*16(ZZ,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPA 
COMPLEX* 16 ZZ, DFDZ , I,RV 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC, KAPPA /FOURTH/ PI 
C 

I = (0 . DO , 1 .DO) 

DELTAT = 1 . DO/DFLOAT ( ISTEPS ) 

TT = TE * 2. DO* PI 



C 
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RV = ZZ - l.DO / DCONJG(ZZ) 

C 

DFDZ -DSIN (TT ) * ( 1 . DO-1 . DO/ (ZZ*ZZ ) ) - I*KAPPA/RV 
C 

Q = DCONJG(DFDZ) 

C 

RETURN 

END 

C 

C 

(2 ********************************************************************* 
C * SUBROUTINE OUTPUT * 

Q ********************************************* ************************ 

C 

C THIS SUBROUTINE IS USED TO SAVE THE FINAL VALUES OF Z,TT, AND NUMCYC 
C FOR POSSIBLE CONTINUATION OF THE PLOTS. 

C THE VALUES ARE SAVED IN FILE 'FILE FT25F001 1 . 

C 

SUBROUTINE OUTPUT (Z) 

IMPLICIT REAL*8 (A-H, O-Z) 

COMPLEX*16 Z(14409) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 
C 

DELTAT = 1 . DO/FLOAT ( I STEPS ) 

NN = NUMPTS - 4 
TTF = DFLOAT(NN)* DELTAT + TTI 
NUMCYC = NN / ISTEPS 
C 

REWIND 25 
DO 10 N = 1,4 
NN = NN + 1 
10 WRITE (25 , *) Z (NN ) 

WRITE (25 , *) TTF 
'WRITE (25 , *) NUMCYC 
C 
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RETURN 

END 

C • 

C 

C * SUBROUTINE FORCE * 

c 

C THIS SUBROUTINE IS USED TO CALCULATE THE FORCE COEFFICIENTS CD, AND 
C CL FOR THE CYLINDER BASED ON THE LOCATIONS OF THE VORTEX AND A 
C FINITE DIFFERENCING SCHEME. 

C 

SUBROUTINE FORCE (Z , CD , CL) 

IMPLICIT REALMS (A-H, O-Z) 

REAL*4 CD(14409) ,CL(14409) 

REALMS KC, KAPPA 
C0MPLEX*8 COEF 
COMPLEX*16 Z(14409) , I 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC, KAPPA 
1 /FOURTH/ PI 

C 

I = (0 . DO , 1 . DO ) 

DELTAT = 1 .DO/DFLOAT(ISTEPS) 

C 

DO 10 N = 2, NUMPTS 
C 

TT = (DFLOAT(N-l) * DELTAT + TTI ) * 2. DO * PI 
C 

COEF = (2 .DO*PI*PI/KC)*DCOS (TT) 

1 + ((PI*KAPPA)/(KC*DELTAT))*I* 

2 ( (Z(N)-l .DO/DCONJG(Z(N) ) ) - 

3 (Z(N-l) - l.DO/DCONJG(Z(N-l)))) 

COEF = CONJG(COEF) 

CD(N) = REAL (COEF) 

10 CL (N) = I MAG (COEF) 
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c 

RETURN 

END 

C 

C 

C ****** ***K**KK**XK7SK*K*7t*k&Xx&***KKk***k****-k***7<:***? W^X^7C7!7tic7C^^7C?: 

C * SUBROUTINE GRAPH * 

AAXX^XXAX^XAXXXAAXXXAXA^AXXAXXXAXAAXAAAXAXA^ 

C 

SUBROUTINE GRAPH (Z,QVECT,DISP,CD,CL) 

IMPLICIT REALMS (A-H, 0-Z) 

REAL*4 ZR( 14409) ,ZI(14409) ,X(361) ,Y(361) ,QMAG(14409) 

REAL*4 DISP(14409) ,CD(14409) ,CL(14409) ,T(14409) 

REAL*4 U( 14409) ,V( 14409 ) ,GARRAY(1) 

COMPLEX*16 Z( 14409) ,QVECT (14409) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /FOURTH/ PI 
C 

C DEFINE CYLINDER OF RADIUS 1.0 
C 

DO 10 N = 1,361 

RAD = FLOAT (N- 1 ) * PI /180.0 

X(N) = COS (RAD) 

10 Y(N) = SIN (RAD) 

C 

C DEFINE REAL SINGLE PRECISION NUMBERS FOR PLOT ROUTINES 
C 

DELTAT = 1 .D0/DFLOAT(ISTEPS) 

C 

DO 20 N =1, NUMPTS 
ZR(N) =DREAL(Z (N) ) 

ZI(N) =DIMAG(Z(N) ) 

QMAG(N) = ABS (QVECT(N) ) 

U(N) = DREAL(QVECT (N) ) 

20 V(N) = DIMAG(QVECT (N) ) 

C 
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C GENERATE PLOTS 



C 



C 

c 

c 

c 

c 

c 

c 



c 



c 



c 



CALL COMPRS 
CALL TEK618 
CALL BLOWUP ( .3365655) 

CALL BLOWUP (.55) 

CALL BLOWUP (2.0) 

CALL NOBRDR 

VORTEX TRACK PLOT 

CALL GRACE (0.0) 

CALL NOCHEK 

CALL PAGE (11.0,3.5) 

CALL AREA2D ( 10 . , 5.71425) 

CALL YNAME ( 1 Y$ ' ,100) 

CALL XNAME ( 1 X$ 1 ,100) 

CALL YAXANG(O.O) 

CALL INTAXS 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF (-7, 1,7, -4, 1,4) 

CALL CURVE(X,Y, 361,0) 

GARRAY ( 1 ) =0.08 

CALL SHADE (X,Y,361,45, GARRAY ,1,0,0) 



CALL CURVE (ZR,ZI,NUMPTS, 20000) 
C CALL CURVE (ZR, ZI ,NUMPTS , 0 ) 

CALL ENDPL(O) 

C 

C MAKE CONTINUATION SELECTIONS 
C 



WRITE (6,*) 1 CONTINUATION MEMO' 
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WRITE (6,*‘) 

WRITE (6 , *) 

WRITE (6,*) 'TERMINATION OF PROGRAM O' 

WRITE (6,*) 'RERUN PROGRAM 1' 

WRITE (6 , *) ' CONTINUE WITH NORMAL SEQUENCE OF PLOTS 2' 



WRITE (6 , *) 

WRITE (6,*) 'PLEASE ENTER YOUR SELECTION: (0,1,2)' 

READ(5,*)L1 
IF(Ll.LT.l)GOTO 100 
IF(Ll.EQ.l)GOTO 100 
CONTINUE 
C 

PHSHFT = TTI - DFLOAT(INT(TTI) ) 

C 

DO 40 N = 1 , NUMPTS 

40 T(N) = (DFLOAT(N-l) * DELTAT + TTI - PHSHFT) 

C 

WRITE (6,*) 'DO YOU WANT A FORCE PLOT? 1-YES, 0-N0 1 
READ (5,*) LI 
IF (Ll.NE.l)GOTO 50 
C 

C FORCE VS TIME PLOTS 
C 

WRITE (6,*) 'ENTER VALUE OF THE LEFT HAND SIDE OF THE TIME SCALE: 
READ (5,*) ILHS 

WRITE (6,*) 'ENTER VALUE OF THE RIGHT HAND SIDE OF THE TIME SCALE 
READ (5 , *) IRHS 
C 

CALL AREA2D(10 . 0 , 7.5) 

CALL YNAME( 'COEF$' ,100) 

CALL XNAME ( ' T*$ ' ,100) 

CALL YAXANG (0.0) 

CALL INTAXS 
C 

CALL GRAF ( ILHS , . 5 , IRHS ,-4,1. ,4.0) 
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c 

CALL CURVE (T, CD, NUMPTS,0)' 

CALL DASH 

CALL CURVE (T, CL, NUMPTS,0) 

CALL RESET ('DASH') 

CALL ENDPL(O) 

C 

50 WRITE (6 ,*) ' DO YOU WANT A Q-MAG VS DISP PLOT? 1-YES, 0-N0 
READ(5 , *) LI 
IF (Ll.NE.l)GOTO 60 
C 

C MAGNITUDE OF Q VS . DISPLACEMENT PLOT 
C 

CALL AREA2D(10 .0 , 7.5) 

CALL YNAME ( ' Q$ 1 ,100) 

CALL XNAME ( ‘ S$ 1 ,100) 

CALL YAKANG(O.O) 

CALL INTAXS 
C 

CALL GRACE (2.) 

CALL GRAF( 0,10, 160, 0,1. ,10.0) 

C 

CALL CURVE (DISP, QMAG , NUMPTS , 0 ) 

CALL ENDPL(O) 

C 

60 WRITE (6 ,*) 1 DO YOU WANT A Q-MAG VS TIME PLOT? 1-YES, O-NO 
READ (5,*) LI 
IF (Ll.NE.l)GOTO 70 
C 

C MAGNITUDE OF Q VS . TIME PLOT 
C 

CALL AREA2D(10 .0 , 7.5) 

CALL YNAME ( 1 Q$ 1 ,100) 

CALL XNAME ( 1 T* $ ' ,100) 

CALL YAXANG(O.O) 
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CALL INTAXS 



C 

CALL GRAF (0,1,10,0,1. ,10.0) 

C 

CALL CURVE (T , QMAG , NUMPTS , 0 ) 

CALL ENDPL(O) 

C 

70 WRITE (6,*) 'DO YOU WANT A U-VEL VS X-POS PLOT? 1-YES, 0-N0 
READ ( 5 , * ) LI 
IF (Ll.NE.l)GOTO 80 
C 

C U-VELOCITY POSITION PLOTS 
C 

CALL AREA2D(10 . 0 , 7.5) 

CALL YNAME ( 1 U-VELOCITY$ 1 , 100) 

CALL XNAME ( 1 X-POSITION$ 1 , 100 ) 

CALL YAXANG (0.0) 

CALL INTAXS 
C 

CALL GRAF(-7,1,7,-10,2. ,10.0) 

C 

CALL CURVE ( ZR , U , NUMPTS , 0 ) 

CALL ENDPL(O) 

C 

80 WRITE (6 , *) 'DO YOU WANT A V-VEL VS Y-POS PLOT? 1-YES, 0-N0 
READ (5,*) LI 
IF (Ll.NE.l)GOTO 90 
C 

C V-VELOCITY POSITION PLOTS 
C 

CALL AREA2D( 10 . 0 , 7.5) 

CALL YNAME ( 1 V-VELOCITY$ 1 ,100) 

CALL XNAME ( 'Y-POSITION$' ,100) 

CALL YAXANG (0.0) 

CALL INTAXS 
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CALL GRAF (-4, .5, 4, -7,1. ,7.0) 

C 

CALL CURVE ( ZI ,V,NUMPTS,0) 

CALL EMDPL(O) 

C 

90 WRITE (6 ,*) 1 DO YOU WANT A U-VEL VS V-VEL PLOT? 1-YES, O-NO' 
READ ( 5 , *) LI 
IF (LI .NE . 1)G0T0 100 
C 

C U VS V VELOCITY PLOTS 
C 

CALL AREA2D (7 . 5 , 7.5) 

CALL YNAME ( 1 V-VELOCITY$ 1 , 100 ) 

CALL XNAME( 1 U-VELOCITYS 1 , 100) 

CALL YAXANG(O.O) 

CALL INTAXS 
C 

CALL GRAF(-10,2,10,-10,2. ,10.0) 

C 

CALL CURVE (U ,V ,NUMPTS , 0) 

CALL ENDPL(O) 

C 

100 CONTINUE 
C 

RETURN 

END 
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APPENDIX B 

RANKINE VORTEX COMPUTER MODEL 



c 

C THIS FUNCTION SUBROUTINE IS WRITTEN FOR USE WITH PROGRAM 'VORTEX 
C FORTRAN 1 . IT CONTAINS VELOCITY CALCULATIONS FOR A RANKINE 
C VORTEX MODEL. IT MAY BE INSERTED IN PLACE OF THE IDEAL 
C VORTEX MODEL (FUNCTION SUBROUITNE Q IN 'VORTEX FORTRAN') 

C 

C 

c - - - 

C | FUNCTION SUBPROGRAM Q (RANKINE MODEL) 

C 

c 

C THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF Q 
C 

COMPLEX FUNCTION Q*16(ZZ,TE) 

IMPLICIT REAL *8 (A-H.O-Z) 

REAL *8 KC, KAPPA 

COMPLEX*16 ZZ, DFDZ , I, RV, VINDUC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC, KAPPA /FOURTH/ PI 
C 

RCORE = 0.1 DO 
I = (0 .DO , 1 . DO) 

DELTAT = 1 .DO/DFLOAT(ISTEPS) 

TT = TE * 2. DO * PI 
C 

RV = ZZ - 1 . DO / DCONJG(ZZ) 

CHECK = ABS (RV) 

C 

IF ( RCORE. GT. CHECK) GOTO 10 
C 

VINDUC = -1 . DO / RV 
GOTO 20 
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c 

10 VINDUC = -RV / RCORE**2 
C 

20 DFDZ =DSIN(TT)*(1.D0-1.D0/(ZZ*ZZ))+KAPPA*I*VINDUC 

c 

Q = DCONJG(DFDZ) 

C 

RETURN 

END 
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APPENDIX C 

GAUSSIAN VORTEX COMPUTER MODEL 



c 

C THIS FUNCTION SUBROUTINE IS WRITTEN FOR USE WITH PROGRAM 'VORTEX 
C FORTRAN'. IT CONTAINS VELOCITY CALCULATIONS FOR A GAUSSIAN TIME 
C DEPENDENT VORTEX MODEL. IT MAY BE INSERTED IN PLACE OF THE IDEAL 
C VORTEX MODEL (FUNCTION SUBROUITNE Q IN 'VORTEX FORTRAN') 

C 

C 

C 

C | FUNCTION SUBPROGRAM Q (GAUSSIAN MODEL) 

C - 

c 

C THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF Q 
C 

COMPLEX FUNCTION Q*16(ZZ,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL *8 KC, KAPPA 
COMPLEX*16 ZZ, DFDZ , I,RV 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC, KAPPA /FOURTH/ PI 
C 

I = (0 .DO , 1 . DO ) 

DELTAT = 1 . DO/DFLOAT ( ISTEPS ) 

TT = TE * 2. DO* PI 

TVINV = 1 . DO /(TE + .497440685D0) 

C 

RV = ZZ - 1 .DO / DCONJG(ZZ) 

RVMAG = ABS(RV) 

ARG = -62.5D0*TVINV*RVMAG**2 
IF (ARG. LE. 50. DO) GOTO 10 
ARG = 50. DO 
C 

10 VMOD = 1 . DO - DEXP (ARG) 
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c 

DFDZ =DSIN(TT)*(1.DO-1.DO/(ZZ*ZZ))-I*(KAPPA/RV)*VMOD 
C 

Q = DCONJG(DFDZ) 

C 

RETURN 

END 
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APPENDIX D 

TWO VORTEX COMPUTER MODEL 






C * * 

C* • TWOVOR * 

c * * 

C * BY LT. W.T. MCCOY * 

C * * 



Q ***************** I*:******** ***************&****&******* *************** 

c 

C DESCRIPTION 

C 

C THIS PROGRAM TRACKS TWO VORTICIES IN COMPLEX SPACE AROUND A UNIT 
C CYLINDER IN OSCILLATING FLOW. INITIAL LOCATION, STARTING TIME, AND 
C THE NUMBER OF CYCLES OF THE OSCILLATING FLOW ARE REQUIRED INPUT. 

C 

C THE PROGRAM IS DESIGNED TO RUN ON DISSPLA BY USING THE FOLLOWING 
C DISSPLA EXEC COMMAND :' DISSPLA TWOVOR'. IT IS CURRENTLY 
C DIMENSIONED TO RUN A MAX OF 6 CYCLES AT 1440 STEPS PER CYCLE. 

C 2.0M OF VIRTUAL MEMORY IS REQUIRED TO RUN THE PROGRAM AT MAXIMUM 
C CAPACITY. THE PROGRAM IS INTERACTIVE FOR INPUT OPTIONS. 

C 

C ALL CONSTANTS ARE SET IN THE MAIN PROGRAM. PREVIOUS INPUT OPTIONS 
C ARE READ FROM FILE 'FILE FT35F001 ' . FAILURE TO GENERATE THIS FILE 
C PRIOR TO INITIAL RUNNING RESULTS IN A ERROR MESSAGE AND TERMINATION 
C OF THE RUN. PRIOR TO RUNNING, FILE 'FILE FT35F001' SHOULD BE 
C INITIATED CONTAINING 8 COMPLEX NUMBERS, 1 REAL NUMBER, AND 1 INTEGER. 
C WITH THIS INPUT THE PROGRAM DOES NOT NEED TO BE TO BE RECOMPILED 
C FOR EACH RUN. 

C 

C THIS PROGRAM IS WRITTEN IN DOUBLE PRECISION AND EMPLOYS THE MILNE 
C PREDICTOR-CORRECTOR CONVECTION SCHEME, WITH INITIAL POINTS GENER- 
C ATED BY A 4TH ORDER RUNGE-KUTTA SCHEME. 
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c 

C OUTPUT OF THE PROGRAM INCLUDES THE LAST FOUR VALUES OF THE LOCATION 
C OF THE VORTICES, FINISH TIME, AND THE NUMBER OF CYCLES. THESE VALUES 
C ARE WRITTEN IN FILE 'FILE FT35F001' AND ALLOW FOR OBSERVATIONS BEYOND 
C THE MAX LIMIT OF 6 CYCLES BY CONTINUING WITH ADDITIONAL RUNS. 

C 

C OUTPUT FOR THE VORTEX PATHS ARE MADE GRAPHICALLY. PLOTS MAY BE MADE 
C ON A TEK618 OR SENT TO A METAFILE (PRE-EDIT SUBROUTINE 'GRAPH'). 

C 

c 

c 

c 

C ALPHABETICAL LISTING OF VARIABLES 

C 

C VARIABLE (SUBROUTINE USED) DESCRIPTION. 

C 

c 

c 

C DELTAT (STEP, Q, OUTPUT, FORCE, GRAPH) INCREMENTAL STEP 

C IN NON-DIMENSIONAL TIME (INVERSE OF STEPS PER 

C CYCLE). 

C 

C DFDZ (QA,QB) COMPLEX PARTIAL DERIVATIVE OF THE COMPLEX 

C FUNCTION F WITH TO RESPECT Z. 

C 

C I (QA,QB, FORCE) CONSTANT COMPLEX NUMBER = (0.0, 1.0) 

C 

C ISTEPS (ALL)NUMBER OF INCREMENTAL STEPS PER CYCLE TO BE 

C TAKEN FOR ALL ITERATIONS. 

C 

C Kl-4 (RUNGE) RUNGE-KUTTA CONSTANTS 

C 

C KAPPAA (ALL) STRENGTH OF THE VORTEX (GAMMA/ (2*PI*UMAX) ) 

C 

C KAPPAB (ALL) STRENGTH OF THE VORTEX (GAMMA/ (2*PI*UMAX) ) 
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c 

C KC (ALL) KEULLEGAN-CARPENTER NUMBER FOR THE OSCILLATING 

C FLOW. 

C 

C LI (INPUT) LOGICAL OPERATOR. 

C 

C NUMCYC (INPUT) TOTAL NUMBER OF CYCLES TO BE EVALUATED. 

C 

C NUMPTS (ALL) TOTAL NUMBER OF POINTS TO BE DETERMINED. 

C (NUMCYC * ISTEPS) 

C 

C PI (ALL)CONSTANT = 3.141592654 

C 

C QA (QA) FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 

C VELOCITY OF THE VORTEX. 

C 

C Q3 (QB)FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 

C VELOCITY OF THE VORTEX. 

C 

C QVECTA(N) (STEP, GRAPH) COMPLEX ARRAY CONTAINING THE VELOCITY 

C OF THE VORTEX AT ANY POINT Z(N). 

C 

C QVECTA(N) (STEP, GRAPH) COMPLEX ARRAY CONTAINING THE VELOCITY 

C OF THE VORTEX AT ANY POINT Z(N). 

C 

C TE (RUNGE ,QA,QB)NON-DIMENSIONAL ENVIRONMENTAL FLOW TIME 

C 

C TT ( STEP , QA , QB ) NON- DIMENSIONAL TIME CHANGED TO 

C RADIANS FOR COS/SIN EVALUATIONS (0-2*PI FOR A CYCLE) 

C 

c TTF (OUTPUT) NON- DIMENSIONAL FINISH TIME OF THE RUNS. 

C 

C TTI (ALL) INITIAL NON-DIMENSIONAL STARTING TIME FOR THE 

C OSCILLATING FLOW. 

C 
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C ZA (ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 

C 

C ZB (ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 

C 

C ZAI (GRAPH) IMAGINARY (Y) LOCATIONOF THE VORTEX. 

C 

C ZBI (GRAPH) IMAGINARY(Y) LOCATIONOF THE VORTEX. 

C 

C ZP (QA,QB) COMPLEX PREDICTED VALUE FOR Z(N) 

C 

C ZAR ( GRAPH) REAL (X) LOCATION OF THE VORTEX. 

C 

C ZBR ( GRAPH) REAL (X) LOCATION OF THE VORTEX. 

C 

C ZZ (QA,QB) DUMMY ARGUMENT FOR THE FUNCTION SUBROUTINE QA 

C AND QB WHICH CONTAINS THE PASSED VALUE OF Z(N) . 

C (RUNGE) DUMMY ARGUMENT CALL FUNCTION SUBROUITNE QA 

CAND QB 

C 

c 

c 

c 



C********************************************************************* 

C * MAIN PROGRAM * 

Q-k'k'k-k'k-k 'k-k’k'k ■k'k'k’k'k’k 

c 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC , KAPPAA , KAPPAB 
COMPLEX*16 ZA(8649) ,QVECTA(8649) 

COMPLEX*16 ZB (8649) ,QVECTB(8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC , KAPPAA, 
l KAPPAB/ FOURTH/ PI 

C 

C SET CONSTANTS 
C 
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KAPPAA = 0.5D0 
KC = 10. DO 
PI = 3 . 141592654D0 
ISTEPS = 1440 
C 

10 CALL INPUT (ZA,Z3) 

CALL STEP (ZA,QVECTA, ZB ( QVECTB) 

CALL OUTPUT (ZA,Z3) 

CALL GRAPH (ZA , QVECTA , Z3 , QVECTB ) 

C 

WRITE (6,*) 'DO YOU WISH TO RE-RUN THE PROGRAM? 1-YES, O-NO 1 
READ (5,*) LI 
IF (Ll.NE.l) GOTO 20 

WRITE (6,*) 'EXECUTION CONTINUES 1 

GOTO 10 
20 CONTINUE 
C 

CALL DONEPL 
STOP 
END 
C 

c 

C * SUBROUTINE INPUT * 

C 

C THIS SUBROUTINE INTERACTIVELY INPUTS THE VALUES FOR THE INITIAL 
C LOCATION OF THE VORTEX, STARTING TIME, ITERATION STEP SIZE, AND 
C THE NUMBER OF CYCLES TO BE EXAMINED. 

C 

SUBROUTINE INPUT (ZA,Z3) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL *8 KC, KAPPAA, KAPPAB 
COMPLEXES ZA(8649) ,ZB(8649) 

COMMON /FIRST/ TTI, ISTEPS /SECOND/ NUMPTS/THIRD/KC , KAPPAA , KAPPAB 
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c 

C READ LAST 4 POINTS OF EACH ARRAY FOR CONTINUATION OPTIONS 
C 

REWIND 35 
DO 10 N = 1,4 
10 READ (35 , *) ZA(N),ZB(N) 

READ (35,*) TTI 
RE AD (3 5,*) NUMCYC 
C 

C CONTINUATION OPTIONS 
C 

WRITE (6,*) 'DO YOU WISH TO CONTINUE THIS RUN FORM THE LAST? 
# O-NO' 

READ (5,*) LI 
IF (L1.NE.0) GOTO 20 
C 

WRITE (6,*) 'INPUT THE VALUE OF X : ZA(1)=(X,X)' 

READ (5,*) X 
ZA(1) = DCMPLX(X,X) 

C 

WRITE ( 6, *)•' INPUT THE VALUE OF EPSILON RADIUS:' 

READ ( 5 , *) ER 

X = ER * X 

ZB ( 1 ) = DCMPLX(X, -X) 

C 

WRITE ( 6 , *) 'INPUT THE VALUE OF EPSILON KAPPA:' 

READ (5,*) EK 
KAPPAB = KAPPAA * EK 
C 

TTI = 0.5D0 
C 

C CALL SUBROUTINE RUNGE TO GENERATE FIRST 4 POINTS FOR EACH ARRAY 
C 

CALL RUNGE (ZA,ZB) 

C 



1-YES 
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C INPUT THE NUMBER OF CYCLES TO BE ANALYZED 
C 

20 WRITE (6,*) 'THE LAST VALUE FOR NUMCYC =' / NUMCYC 

WRITE(6 , *) 1 DO YOU WISH TO CHANGE IT? 1-YES O-NO 1 
READ (5,*) LI 
IF(Ll.NE.l) GOTO 30 

WRITE (6 ,*) 'INPUT THE NEW VALUE OF NUMCYC: 1 
READ (5,*) NUMCYC 
30 CONTINUE 
C 

NUMPTS = NUMCYC * ISTEPS + 4 
C 

RETURN 

END 

C 

C 

C * SUBROUTINE RUNGE * 

Q A****** A***** Atc* ** *** *** kirk A************** 

c 

C THIS SUBROUTINE USES A 4TH ORDER RUNGE-KUTTA METHOD TO GENERATE THE 
C FIRST FOUR VALUES OF Z. 

C 

SUBROUTINE RUNGE (ZA, ZB) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX*16 ZA(8649) ,ZB(8649) ,QA,K1 ,K2 ,K3 ,K4,ZP1 ,ZZ 
COMMON /FIRST/ TTI, ISTEPS 
C 

DELTAT = 1 .DO/DFLOAT( ISTEPS) 

C 

C RUNGE-KUTTA EVALUATION FOR FIRST FOUR POINTS 
C 

DO 5 N = 1,3 
’ ZZ = ZA(N) 

TE1 = (DFLOAT(N-l) * DELTAT + TTI) 
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K1 = DLETAT * QA(ZZ,ZB(N) ,TE1) 

22 = ZA(N)+0 . 5D0*K1 

TE2 = TE1+ 0 . 5D0*DELTAT 

K2 = DELTAT * QA(ZZ,ZB(N) ,TE2) 

ZZ = ZA(N)+0.5D0*K2 

K3 = DELTAT * QA(2Z ; 2B(N) > TE2) 

ZZ = ZA (N) + K3 

TE4 = TE1 + DELTAT 

K4 = DELTAT * QA(ZZ,ZB(N) / TE4) 

C 

ZA(N+1) = ZA(N) + (K1-2.D0*K2+2.D0*K3+K4)/6.D0 
C 

ZZ = ZB(N) 

K1 = DLETAT * QB(ZA(N) ,ZZ,TE1) 

ZZ = ZB(N)+0 . 5D0*K1 

K2 = DELTAT * QB(ZA(N) ,ZZ ,TE2) 

ZZ = ZB(N)+0.5D0*K2 

K3 = DELTAT * QB(ZA(N),ZZ, TE2 ) 

ZZ = ZB (N) + K3 

TE4 = TE1 + DELTAT 

K4 = DELTAT * QB(ZA(N) ,ZZ,TE4) 

C 

5 ZB (N+l ) = ZB(N) + (K1-2.D0*K2+2.D0*K3+K4)/6.D0 
C 
C 

RETURN 

END 

C 

C 

C * SUBROUTINE STEP * 

Q ***** ****** K * 7* 

C 

C THIS SUBROUTINE PROGRESSIVELY CALCULATES THE VALUES OF ZA(N) AND 
C ZB(N) BASED ON PREDICTOR-CORRECTOR CONVECTION SCHEME. 
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c 

SUBROUTINE STEP ( ZA , QVE CTA , ZB , QVE CTB ) 

IMPLICIT REAL* 8 (A-H,0-Z) 

REAL*8 KC , KAPPAA , KAPPAB 

COMPLEX* 16 ZA(8649) ( QVECTA(8649) ,ZP,QA,QB 
COMPLEX*16 ZB(8649) ,QVECTB(8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS/THIRD/ KC, KAPPAA, KAPPAB 
C 

DELTAT = 1 . DO/DFLOAT ( I STEPS ) 

C 

C MAKE VELOCITY EVALUATIONS FOR POINTS 1 THROUGH 4 
C 

DO 7 N = 1,4 

TE = (DFLOAT (N-l ) * DELTAT + TTI ) 

QVECTA(N) = QA(ZA(N) ,ZB(N) ,TE) 

7 QVECTB(N) = QB(ZA(N) ,ZB(N) ,TE) 

C 

DO 10 N = 5, NUMPTS 

NM1 = N-l 

NM2 = N-2 

NM3 = N-3 

NM4 = N-4 

TE = (DFLOAT(N-l)*DELTAT+TTI) 

C 

C CALCULATE PREDICTED VALUE OF ZA(N) BASED ON ZB(N-l) 

C 

ZP = ZA(NM1) + (55.D0*QVECTA(NM1) - 59 . D0*QVECTA (NM2 ) 

1 +37 . D0*QVECTA(NM3 ) -9 . D0*QVECTA(NM4 ) ) 

2 *DELTAT* 2. DO * KC / 24. DO 
C 

C CALCULATE 1ST ESTIMATE OF ZA(N) BASED ON ZB (N-l) 

C 

ZA(N) = ZA(NMl) + (9 . D0*QA (ZP , ZB (NM1 ) ,TE) + 19 . D0*QVECTA(NM1 ) 

1 -5.D0*QVECTA(NM2)+QVECTA(NM3) )*DELTAT * 2. DO * KC / 24. DO 

C 
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C CALCULATE PREDICTED VALUE OF ZB(N) BASED ON ESTIMATED ZA(N) 

C 

ZP = ZB (NM1 ) + ( 55 . D0*QVECTB (NM1 ) -59 . D0*QVECTB (NM2 ) 

1 +37 . D0*QVECTB (NM3 ) - 9 . D0*QVECTB(NM4 ) ) 

2 *DELTAT* 2. DO * KC / 24. DO 
C 

C CALCULATE 1ST ESTIMATE OF ZB(N) BASED ON ESTIMATED ZA(N) 

C 

ZB (N) = ZB(NMl) + (9.D0*QB(ZA(N) ,ZP,TE) + 19 .D0*QVECT3 (NM1 ) 

1 -5 .D0*QVECTB (NM2 )+QVECTB(NM3 ) )*DELTAT * 2. DO * KC / 24. DO 

C 

C CALCULATE PREDICTED VALUE OF ZA(N) BASED ON ESTIMATED ZB(N) 

C 

ZP = ZA(NM1)+(55.D0*QVECTA(NM1) -59.DO*QVECTA(NM2) 

1 +37 . D0*QVECTA(NM3) -9 . D0*QVECTA(NM4) ) 

2 *DELTAT* 2. DO * KC / 24. DO 
C 

C CALCULATE ACTUAL ZA(N) BASED ON ESTIMATED ZB(N) 

C 

ZA(N) = ZA(NMl) + (9.D0*QA(ZP,ZB(N) ,TE) + 19 .D0*QVECTA(NM1 ) 

1 - 5 . DO*QVECTA (NM2 ) +QVECTA (NM3 ) )*DELTAT * 2. DO * KC / 24. DO 

C 

C CALCULATE PREDICTED VALUE OF ZB(N) BASED ON ACTUAL ZA(N) 

C 

ZP = ZB(NMl)+(55 . D0*QVECTB (NM1 ) -59 .D0*QVECTB (NM2) 

1 +37 . D0*QVECTB (NM3 ) -9 . DO*QVECTB (NM4 ) ) 

2 *DELTAT* 2. DO * KC / 24. DO 
C 

C CALCULATE ACTUAL ZB(N) BASED ON ACTUAL ZA(N) 

C 

ZB (N) = ZB(NMl) + (9 .D0*QB(ZA(N) ,ZP,TE) + 19 . DO*QVECTB(NMl ) 

1 -5.D0*QVECTB(NM2)+QVECTB(NM3) )*DELTAT * 2. DO * KC / 24. DO 

C 

C EVALUATE VELOCITIES AT ZA(N) AND ZB(N) 

C 
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QVECTA (N) = QA(ZA(N) ,ZB(N) ,TE) 

10 QVECTB(N) = QB(ZA(N) ,ZB(N) ,TE) 

C 

RETURN 

END 

C 

C 

C 

C | FUNCTION SUBPROGRAM QA 

C - 

c 

C THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF QA 
C 

COMPLEX FUNCTION QA*16(ZA,ZB,TE) 

IMPLICIT REAL*3 (A-H,0-Z) 

REAL*8 KC, KAPPAA , KAPPAB 
COMPLEX* 16 ZA,ZB, DFDZ , I,AA,AB,AC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC , KAPPAA , KAPPAB /FOURTH/ PI 
C 

I = (0 .DO , 1 .DO) 

TT = TE * 2. DO* PI 
C 

AA = 1 .D0/(ZA-1 .DO/ (DCONJG(ZA) ) ) 

AB = 1 .DO/ (ZA-ZB) 

AC = 1 .DO/ (ZA-1 .DO/ (DCONJG(ZB) ) ) 

C 

DFDZ =DSIN(TT)*(1.D0-1.D0/(ZA*ZA))-I*KAPPAA*AA-I*KA?PA3*(AE-AC) 
C 

QA = DCONJG(DFDZ) 

C 

RETURN 

END 

C 



C 
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FUNCTION SUBPROGRAM QB 



C - 
C I 
C - 

C 

C THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF QB 
C 

COMPLEX FUNCTION QB*16 (ZA , 2B , TE ) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPAA , KAPPAB 
COMPLEX*16 ZA,ZB , DFDZ , I,BA ; BB,BC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC , KAPPAA , KAPPAB /FOURTH/ PI 
C 

I = (0 .DO , 1 .DO) 

TT = TE * 2. DO* PI 
C 

BA = 1 .DO/ (ZB-ZA) 

BB = 1.D0/ (ZB-1.D0/ (DCONJG(ZA) ) ) 

BC = 1 .D0/(ZB-1 .DO/ (DCONJG(ZB) ) ) 

C 

DFDZ =DSIN(TT)*(1 .DO-1 .DO/ (ZB*ZB) )+I*KAPPAA*(BA-BB)+I*KAPPAB*BC 
C 

QB = DCONJG(DFDZ) 

C 

RETURN 

END 

C 

C 

Q ********************************************************************* 

C * SUBROUTINE OUTPUT * 

Q ********************************************************************* 

c 

C THIS SUBROUTINE IS USED TO SAVE THE FINAL VALUES OF Z,TT, AND NUMCYC 
C FOR POSSIBLE CONTINUATION OF THE PLOTS. 

C THE VALUES ARE SAVED IN FILE 'FILE FT35F001 1 . 

C 
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SUBROUTINE OUTPUT(ZA ; ZB) 

IMPLICIT REAL*8 (A-H, O-Z) 

C0MPLEX*16 ZA( 8649), ZB (8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 
C 

DELTAT = 1. DO/FLOAT (ISTEPS) 

NN = NUMPTS - 4 
TTF = DFLOAT(NN)* DELTAT + TTI 
NUMCYC = NN / ISTEPS 
C 

REWIND 35 
DO 10 N = 1,4 
NN = NN + 1 

10 WRITE (35,*) ZA(NN) ,Z3(NN) 

WRITE (35,*) TTF 
WRITE (3 5,*) NUMCYC 
C 

RETURN 

END 

C 

C 

c 

c 

Q X****X****XX*£ ******************* X****************** ^-fr****** *********: 

C * SUBROUTINE GRAPH - . * 

c 

SUBROUTINE GRAPH ( ZA , QVE CT A , ZB , QVE CTB ) 

IMPLICIT REAL*8 (A-H, O-Z) 

REAL*4 ZRA(8649) ,ZIA(8649) ,X(361) ,Y(361) , GARRAY ( 1 ) 

REAL*4 ZRB(8649) ,ZIB(8649) 

COMPLEX* 16 ZA(8649) ,QVECTA(8649) 

COMPLEX* 16 ZB(8649) ,QVECTB(8649) 

COMMON /FIRST/ TTI, ISTEPS /SECOND/ NUMPTS /FOURTH/ PI 
C 
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C DEFINE CYLINDER OF RADIUS 1.0 
C 

DO 10 N = 1,361 

RAD = FLOAT (N-l) * PI /180.0 

X(N) = COS (RAD) 

10 Y(N) = SIN (RAD) 

C 

C DEFINE REAL SINGLE PRECISION NUMBERS FOR PLOT ROUTINES 
C 

DELTAT = 1 .DO/DFLOAT(ISTEPS) 

NPTSM4 = NUMPTS - 4 
C 

DO 20 N =1, NPTSM4 
ZRA(N) =DREAL(ZA(N) ) 

ZIA(N) =DIMAG(ZA(N) ) 

ZRB(M) =DREAL(ZB (N) ) 

20 ZIB(N) =DIMAG(ZB(N) ) 

C 

C GENERATE PLOTS 
C 

C CALL COMPRS 

C CALL BLOWUP (1.0) 

CALL TEK618 
CALL BLOWUP (2.0) 

CALL NOBRDR 
CALL GRACE (0.0) 

C 

C VORTEX TRACK PLOT 

C 

CALL NO CHE K 

CALL PAGE (11.0,8.5) 

CALL AREA2D(8. 0,6.0) 

CALL YNAME ( 1 Y$ 1 ,100) 

CALL XNAME ( 1 X$ 1 ,100) 

CALL YAXANG (0.0) 
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CALL INTAXS 



C 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF (-26,1,2,-2,1,19) 

C 

CALL CURVE(X,Y, 361,0) 

GARRAY ( 1 ) =0.08 

CALL SHADE(X, Y, 361 ,45 , GARRAY, 1,0,0) 
C 

CALL CURVE ( ZRA , ZIA , NPTSM4 , 0 ) 

C 

CALL DASH 

CALL CURVE ( ZRB , ZIB , NPTSM4 , 0 ) 

CALL ENDPL(O) 

C 

RETURN 

END 
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APPENDIX E 
FIGURES 





Figure E.Ia Experimental In-line Force Trace 
Top - Force CoefTicient, Bottom - Flow Velocity 
K c - 10.0 
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Figure F.lb Experimental Transverse Force Trace 
Top - Force Coefficient, Bottom - Flow Velocity 
K = 10.0 

C 
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Figure F.2 Vortex Path in Uniform Flow 
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Figure E.3 Force CoefFicicnts for Z( 1) = ( 1.9, 1.9) and Ti- 0.75 
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Figure E.4 Force Coefficients for Z(1)=(0, 2.2192) and Ti = 0.75 
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Ficure E.5 Force CocfTicients for Z(l)~- (0 t 2.2193) and Ti 0.75 
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figure E.6 Force Coefficients for 2(1) = (0, 1.6SS5) and Ti = 0.25 
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Figure E.7 Vortex Path Z(1) = (0, 1.50S3) 
Ti = 0.0 

(a) Cycles 1-5, (b) Cycle 6, (c) Cycles 7-15 



71 






X 



Figure E. 8 Vortex Path Z(l)= (0, 1.5084) 
Ti = 0.0 

(a) Cycles 1-5, (b) Cycle 6, (c) Cycles 7-15 
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Figure E.9 Vortex Paths for 100 Cycles 
Ti= 0.0 

(a) Z(1) = (0, 1.5083), (b) Z(1)=(0, 1.5084) 
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Figure E.10 Velocity Magnitudes for 10 Cycles 
Ti = 0.0 

(a) Z(1) = (0, 1.50S3), (b) Z( 1)= (0, 1.5034) 
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Figure E.l 1 u-Velocity vs x-Position for 10 Cycles 

Ti = 0.0 

(a) Z(l) = (0, 1.50S3), (b) Z(1) = (0, 1.5084) 
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Figure E.12 v-Velocity vs y-Position for 10 Cycles 

Ti = 0.0 



(a)Z(l)-(0, 1.5083), (b) Z(l)-(0, 1.5084) 
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Figure E.13 u-Velocity vs v-Velocitv for 10 Cycles 

Ti = 0.0 



(a) Z(1) = (0, 1.50S3), (b) Z(1) = (0, 1.5084) 
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Figure E.14 Vortex Path Z(1) = (0, 1.6SS5) 
Ti = 0.25 

(a) Cycles 1-3, (b) Cycle 4, (c) Cycles 5-15 
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Figure E.15 Vortex Path Z(1) = (0, 1.6886) 
Ti = 0.25 

(a) Cycles 1-3, (b) Cycle 4, (c) Cycles 5-15 
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Figure E.16 Velocity Magnitudes for 10 Cycles 
Ti = 0.25 

(a) Z(1) = (0, 1.6885), (b) Z(l)-(0, 1.6886) 
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Figure E. 17 u-Velocity vs x-Position for 10 Cycles 

Ti = 0.25 



(a) Z(l) = (0, 1.6885), (b) Z(1)»(0, 1.6S86) 
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Figure E.18 v-Velocity vs y-Position for 10 Cycles 

Ti = 0.25 



(a) Z( 1) = (0,1.6885), (b) Z(1) = (0, 1.6886) 
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Figure E.19 u-Velocity vs v-Velocity for 10 Cycles 

Ti = 0.25 

(a) Z(1) = (0, 1.6885), (b) Z(l)-(0, 1.6886) 



83 



4 i 
3- 
2 - 
1 - 

(a) >- 0- 

-1 - 
-2 - 
-3 - 

~ V ' i ' i' 1 i ' I ' l 1 l ' l 1 l " I ' 1 ' I ' l ’ l 1 l 

-7-6-5-4-3 -2 -1 0 1 2 3 4 5 6 7 

X 

4 -i 
3- 
2 - 
1 - 

(b) >, o- 

-1 ~ 

- 2 - 
-3- 

-4 ■ i ' i 1 i ' "i ■ i i i i -i - |— ■ — r--> — | i | i | ■ | i | 

-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 

X 






X 

Figure E.20 Vortex Path Z(1) = (0, 2.2192) 
Ti = 0.75 

(a) Cycles 1-5, (b) Cycle 6, (c) Cycles 7-15 
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Figure E.21 Vortex Path Z(l) = (0, 2.2193) 
Ti = 0.75 

(a) Cycles 1-5, (b) Cycle 6, (c) Cycles 7-15 
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Figure E.22 Velocity Magnitude for 10 Cycles 
Ti = 0.75 



(a) Z(1) = (0, 2.2192), (b) Z(1) = (0, 2.2193) 
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Figure E.23 u-Velocity vs x-Position for 10 Cycles 



Ti = 0.75 

(a) Z(1) = (0, 2.2192), (b) Z(1)=(0, 2.2193) 
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Figure E.24 v-Velocity vs y-Position for 10 Cycles 

Ti = 0.75 

(a).Z( 1) = (0, 2.2192), (b) Z(1)=(0, 2.2193) 
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Figure E.25 u-Velocity vs v-Velocity for 10 Cycles 

Ti = 0.25 

(a) Z(l) = (0, 2.2192), (b) Z(1) = (0, 2.2193) 
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Figure E.26 Vortex Path Z(l)=(1.0, 1.0272) 
Ti = 0.0 

(a) Cycles 1-6, (b) Cycle 7, (c) Cycles 8-15 
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Figure E.27 Vortex Path Z(l) = (1.0, 1.0273) 
Ti = 0.0 



(a) Cycles 1-6, (b) Cycle 7, (c) Cycles 8-15 
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Figure E.28 Velocity Magnitudes for 10 Cycles 
Ti = 0.0 

(a) Z(l) = (1.0, 1.0272), (b) Z(l) = (1.0, 1.0273) 
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Figure E.29 u-Velocity vs x-Position for 10 Cycles 

Ti = 0.0 

(a) Z(l) = (1.0, 1.0272), (b) Z(l) = (1.0, 1.0273) 
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Figure E.30 v-Velocity vs y-Position for 10 Cycles 

Ti = 0.0 

(a) Z(l) = (1.0, 1.0272), (b)Z(l) = (1.0, 1.0273) 
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Figure E.31 u-Velocity vs v-Velocity for 10 Cycles 

Ti = 0.0 

(a) Z(l)= (1.0, 1.0272), (b) Z(l)-(1.0, 1.0273) 
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Figure E. 32 Investigation of Phase Shift and Vortex Path Z(l) = (l,l) 
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Figure E.33 Investigation of Phase Shift and Vortex Path Z(l)~ (0,1.50) 

(a) Ti-0.75, (b) Ti = 0.85, (c) Ti=0.92, (d) Ti= 100 
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Figure E.34 Comparison of Vortex Models 
Z( 1) = (0. 1.0273), Ti= 0.0 

(a) Ideal, (b) Rankine, (c) Gaussian Time Dependent 
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Figure E.35 Symmetric Two Vortex Model 
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Figure E.36 Two Vortex Model - Asymmetric Locations 
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Figure E.37 Two Vortex Model - Asymmetric Locations 




102 



Figure E.38 Two Vortex Model * Asymmetric Locations 
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Figure E.39 Two Vortex Model - Asymmetric Strengths 

K a = 0.50, c.. = 0.S5 
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Figure E.40 Two Vortex Model - Asymmetric Strengths 
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Figure E.41 Two Vortex Model - Asymmetric Strengths 
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Figure E.42 Two Vortex Model - Asymmetric Strengths 

K a = 0.50, c.. = 0.975 
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Figure E.43 Two Vortex Model - Asymmetric Strengths & Locations ^ 
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Figure E.44 Two Vortex Model - Cylinder Influenced Z A = (0.8, 0.8), Zg = (0.8, -0.8) 
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